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Abstract 

Swimming microorganisms rely on effective mixing strategies to achieve efficient nutrient influx. Recent 
experiments, probing the mixing capability of unicellular biflagellates, revealed that passive tracer particles 
exhibit anomalous non- Gaussian diffusion when immersed in a dilute suspension of self- motile Chlamy- 
domonas reinhardtii algae. Qualitatively, this observation can be explained by the fact that the algae 
induce a fluid flow that may occasionally accelerate the colloidal tracers to relatively large velocities. A sat- 
isfactory quantitative theory of enhanced mixing in dilute active suspensions, however, is lacking at present. 
In particular, it is unclear how non- Gaussian signatures in the tracers' position distribution are linked to 
the self-propulsion mechanism of a microorganism. Here, we develop a systematic theoretical description 
of anomalous tracer diffusion in active suspensions, based on a simplified tracer- swimmer interaction model 
that captures the typical distance scaling of a microswimmer's flow field. We show that the experimentally 
observed non-Gaussian tails are generic and arise due to a combination of truncated Levy statistics for the 
velocity field and algebraically decaying time correlations in the fluid. Our analytical considerations are 
illustrated through extensive simulations, implemented on graphics processing units to achieve the large 
sample sizes required for analyzing the tails of the tracer distributions. 



Brownian motion presents one of the most beautiful manifestations of the central limit theorem 
in Nature [l[]. Reported as early as 1784 by the Dutch scientist Jan van Ingenhousz [2], the 
seemingly unspectacular random motion of mesoscopic particles in a liquid environment made an 
unforeseeable impact when Perrin's seminal experiments of 1909 [3] yielded convincing evidence for 
the atomistic structure of liquids. This major progress in our understanding of non-living matter - 
which happened long before direct observations of atoms and molecules came within experimental 
reach - would not have been possible without the works of Sutherland |4( and Einstein [5], who 
were able to link the microscopic properties of liquids to a macroscopic observable, namely the 
mean square displacement of a colloidal test particle. 

Caused by many quasi-independent random collisions with surrounding molecules, Brownian 
motion in a passive liquid is quintessentially Gaussian, as predicted by the central limit theorem. 
Remarkably, however, recent experiments by Leptos et al. [6] revealed notable non-Gaussian fea- 
tures in the probability distribution of a tracer particle, when a small concentration of microscale 
swimmers, in their case unicellular biflagellate Chlamydomonas reinhardtii algae, was added to the 
fluid. Understanding this apparent violation of the central limit theorem presents a challenging un- 
solved problem, whose solution promises new insights into the mixing strategies of microorganisms 
[3]. Here, we shall combine extensive analytic and large-scale numerical calculations to elucidate 
the intimate connection between the flow field of an individual microorganism and the anomalous 
(non-Gaussian) diffusion of tracer particles in a dilute swimmer suspension. 

Modern high-speed microscopy techniques resolve the stochastic dynamics of micron-sized tracer 
particles to an ever increasing accuracy [fil, [J]. This opens the exciting possibility of using high- 
precision tracking experiments to probe the statistics of the flow fields created by active swimmers, 
and hence their connection to physical pro perties and evolutionary strategies of microorganisms 
that live in liquid environments 16, 10, 111. Furthermore, a novel class of micromechanical de- 

nn 

vices [12|, [13[ use nonequilibrium fluctuations generated by bacteria as a fundamental ingredient of 
their operation. To explain and exploit the nonequilibrium conditions in active suspensions, it is 
important to fully understand the relation between the experimentally observed features of tracer 
displacements and the characteristics (such as self-propulsion mechanisms) of the algae or bacteria. 

The observations of Leptos et al [(| demonstrate that the time-dependent probability distri- 
butions of tracer displacements in dilute algae suspensions exhibit tails that decay much more 
slowly than would be expected if the tracers obeyed purely Gaussian statistics. At high swimmer 



concentrations, enhanced transport might be expected 
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171 ], as collective behavior emerges from 



swimmer interactions, which can lead to the formation of large-scale vortices and jets. In dilute 
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suspensions, however, where swimmer-swimmer interactions can be neglected, a satisfactory quan- 
titative understanding of the underlying velocity statistics is still lacking. Below, we are going 
to show that the velocity distribution produced by the swimmers takes a tempered (or regular- 
ized) Levy form Q, [jjj], and that the long-time behavior of the tracers' positional probability 
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distribution function can be understood in terms of correlated truncated Levy flights 

If a microswimmer is self-propelled, with no external forces acting, its flow field scales with 
distance as r~ n for an exponent n > 2 [2l|]. We will demonstrate that it is this form of the 
power-law decay that is responsible for the anomalous diffusion of tracer particles [6]. Remarkably, 
qualitatively different behavior can be expected in suspensions of sedimenting swimmers: If gravity 
plays an important role for the swimmer dynamics, the far field flow decays as r _1 and tracers 
will diffuse normally. Finally, our results suggest that, on sufficiently long times scales, anomalous 
tracer diffusion in dilute active suspensions can be viewed as a natural example of a stochastic 
process described by a fractional diffusion equation. 



I. MODEL 



Given an advecting flow ujs[{t,r), generated by a dilute suspension of N self-swimming mi- 
croorganisms, we model the dynamics x(t) of a passive, colloidal tracer particle (radius a) by an 
overdamped Langevin equation of the form 

j t x{t) = u N (t, x{t)) + y/2D^{t). (1) 

The random function £(£) = (£i(£))i=i,2,3 represents uncorrelated Gaussian white noise with 
= an d (0) = SijS(t — £'), describing stochastic collisions between the tracer 

and surrounding liquid molecules. The thermal diffusion coefficient .Do in a fluid of viscosity 77 
is determined by the Stokes-Einstein relation Do = k^T '/ '(Gtvtjcl) . 

If the Reynolds number is very small, the net flow due to a = 1,...,7V active swimmers, 
located at positions X°(t) and moving at velocity V a (t), is, in good approximation, the sum of 
their individual flow fields i£, 

N 

u N (t,x) = J>(*|r*(f)) , r*(f) := (X*(t), V*(t)). (2) 

(7=1 

Since we are interested in physical conditions similar to those in the experiments of Leptos et al [fj], 
our analysis will focus on a dilute suspension of active particles, corresponding to the limit of a small 
volume filling fraction j?« 1. In this case, binary encounters between swimmers are negligible 
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perturbations. Moreover, we can ignore random reorientation of swimmers caused by rotational 
diffusion (due to thermal fluctuations) and search behavior (like chemotaxis or phototaxis), since 
these effects take place on the order of several seconds and thus are not relevant to the tracer 
dynamics. Indeed, in dilute homogeneous solutions, it is irrelevant for the tracer statistics (even on 
longer time scales) whether a tracer experiences two successive scatterings from the same tumbling 
swimmer or from two different non-tumbling swimmers. It is therefore sufficient to assume that 
each swimmer moves ballistically, 

x%t) = xz + tvz. (3) 

For dilute suspensions, the initial swimmer coordinates r cr (0) = (-X"q, Vq) are independent and 
identically distributed random variables with one-particle PDF $i(Tq). More specifically, we 
assume that the distribution of the initial positions Xq is spatially uniform and that the swimmers 
have approximately the same speed V, that is Q>i(Tq) oc 5(V — \ Vq\). 

To complete the definition of the model, we need to specify the flow field u(x\T a (t)) generated by 
a single swimmer. There are various strategies for achieving directed propulsion at the microscale 



221 ]. Small organisms, like algae and bacteria, can swim by moving slender filaments in a manner 
not the same under time reversal. Self-motile colloids, a class of miniature artificial swimmers, 
are powered with interfacial forces induced from the environment Although both of these 

are active particles, microscopic details of their geometry and self-propulsion can lead to different 
velocity fields. This, in turn, affects how a tracer migrates in their flow. In the Stokes regime, if 
external force are absent, self-propelled particles or microorganisms generate velocity fields decaying 
as r~ 2 or faster 
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231 ]. Since we are interested in the general features of mixing by active 
suspensions, and there is no universal description of the flow around an active object, it is helpful 
to consider simplified velocity field models that capture generic features of real microflows. 

We shall focus on two simplified models (see Fig. [TJ that can be interpreted as contributions in 
a general multipole expansion of a flow field. Specifically, we will compare a co-oriented model 
with trivial angular dependence 



\R a \ n + X n 

to a more realistic dipolar (or stresslet) flow field [2 



u(x\T°) = (kV) ( - ) n , n > 1 (4a) 



u(x\T a ) = (kV) 



e 2 



3 ft R) -1 



R . (4b) 



\R a \ 2 + A 2 / 

In Eqs. (jl]), the vector R a {t) :— X a (t) — x connects the swimmer and tracer position at time 
t, R° ' (t) :— R a /\R a \ is the associated unit vector, and (l? :— V° ' /V defines the swimmer's 
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FIG. 1: Flow fields of the model swimmers considered in our analytical and numerical calculations. In both 
plots, the swimmer is oriented upwards and the flow is normalized by the swimming speed. The arrows 
indicate velocity direction and the colors represent magnitude. (Left) Co-oriented model with r~ 2 decay 
and (Right) stroke- averaged flow generated by a dipolar "pusher" . Parameters are those of Fig. [2j 



orientation and swimming direction. The parameter e characterizes the swimmer length scale, k is 
a dimensionless constant that relates the amplitude of the flow field to the swimmer speed, and A 
regularizes the singularity of the flow field at small distances. The co-oriented model (j4aj), due to 
its minimal angular dependence, is useful for pinpointing how the tracer statistics depend on the 
distance scaling of the flow field. For n = 1, the scaling is equivalent to that of an "active" colloid 
or forced swimmer, whereas for n > 2 the scaling resembles that of various natural swimmers not 
subjected to an external force. In particular, the case n — 2 allows us to ascertain the effects of 
the angular dependence of the flow field structure on tracer diffusion, by comparing against the 

atter is commonly considered as a simple stroke-averaged 



more realistic dipolar model (l4bl) . The 
description for natural microswimmers 



As shown in Ref. [24], stroke-averaged models are 
able to capture the most important aspects of the tracer dynamics on time scales longer than the 
swimming stroke of a microorganism. 



II. RESULTS 



We are interested in computing experimentally accessible, statistical properties of the tracer 
particles, such as their velocity PDF, correlation functions, and position PDF. These quantities 
are obtained by averaging suitably defined functions with respect to the TV-swimmer distribution 
— ricr^ifTo)- A detailed description of the averaging procedure and a number of exact 
analytical results are given in the Supplementary Material, below we shall restrict ourselves to 
discussing the main results and their implications. 

We begin by considering the equal-time velocity PDF and velocity autocorrelation function at 
a fixed point in the fluid. Since we are primarily interested in the swimmer contributions, we will 
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focus on the deterministic limit Dq = first. The additive effect of thermal Brownian motion will 
be taken into account later, when we discuss the position statistics of the tracer. Considering a 
suspension of N swimmers, confined by a spherical volume of radius A, the equal-time velocity 
PDF 4>n,a(v) and velocity autocorrelation function Cjv,a(£) of the flow field near the center of the 
container are formally defined by 

<I>n,a( v ) •= (6(v-u N (0,0))), (5a) 
C NtA (t) := (u N (t, 0)u N (0, 0)), (5b) 

where the average ( • ) is taken with respect to the spatially uniform initial distribution of the 
swimmers. For the models (jlj), it is possible to determine c/>n,a and Cjv,a analytically. 



A. Velocity PDF: Slow convergence of the central limit theorem for n > 2 

To elucidate the origin of the unusual velocity statistics in an active suspension, let us consider 
the tracer velocity PDF when there is a single swimmer present, (/>i a(v). The tail of this function 
reflects large tracer velocities generated by close encounters with the swimmer. It is instructive 
to start with the (unphysical) limit A = 0, where the interaction diverges at short distances and 
there is no cutoff for large velocities. In this case, one readily finds from (|5aj) that asymptotically 
0i,a(O,v) oc |v|~( 3 + 3 / n ). This means that the variance of the probability distribution is finite for 
n = 1, but infinite for n > 2. According to Eq. (|2j), the flow field due to N swimmers is the sum 
of N independent and identically distributed random variables. Hence, the central limit theorem 
predicts that, for A = and n = 1, the velocity distribution (j)N,A converges to a Gaussian in the 
large N limit, whereas for A = and n > 2 one expects non-Gaussian behavior due to the infinite 
variance of 0i a- 
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For a real swimmer, the flow field is strongly increasing in the vicinity of the swimmer 
but remains finite due to lubrication effects and nonzero swimmer size. This corresponds in our 
model to a positive value of A. For A > 0, the variance of the one-swimmer PDF (f)i^\ remains finite 
and, formally, the conditions for the central limit theorem are satisfied for all n > 1. However, 
for n > 2 the variance of 0i 5 a remains very large and the convergence to a Gaussian limiting 
distribution is very slow. Our subsequent analysis demonstrates that the velocity PDF is more 
accurately described by a tempered Levy-type distribution. 

These statements are illustrated by Fig. [21 which shows velocity PDFs obtained numerically 
(symbols) and from analytical approximations (solid curves) for the co-oriented model with n — 1 
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(A) and n — 2 (B), and the dipolar model (C) at different swimmer volume fractions tp := 7V(e/A) 3 . 
As evident from Fig. E] A, for n — 1 the velocity PDF indeed converges rapidly to the Gaussian 
shape, in accordance with the central limit theorem. By contrast, for velocity fields decaying as r~ n 
with n > 2, the convergence is surprisingly slow and one observes strongly non-Gaussian features 
at small filling fractions tp. The arrows highlight this regime, which shows a power-law dependence 
of the velocity distribution on the magnitude of the velocity, a signature of a Levy distribution. 

The remarkably slow convergence to the Gaussian central limit theorem prediction can be 
understood quantitatively by considering the characteristic function 

XM) := J d 3 v exp(-iqv) <f>(v) (6) 

of the velocity PDF ([5al) . A detailed analytical calculation (see Supplementary Material) shows 
that the exact result for X(/>(o) can be approximated by 

X,(|g|) = exp{-[(c \q\ 2 + fY' 2 - (7) 

This expression, which reduces to a Gaussian for a — 2, is of the tempered Levy form, and gives 
rise to the following tracer velocity moments 

(| v | 2 ) = 3acfi a ~ 2 , (8a) 
(| v | 4 ) = 15acV~ 4 [2 + a (fi a - 1)] . (8b) 

By studying asymptotic behavior in the small cutoff limit A —> one finds that, for velocity fields 
u decaying as r _n , 

n = 1 a = 2 , n>2 a = 3/n. (9) 

This result confirms that for colloidal- type interactions with n — 1 the velocity PDF is Gaussian, 
whereas for n > 2 deviations from Gaussianity occur in agreement with our numerical results 
of Fig. [2 In the limit [i — 0, Eq. © describes the family of Levy stable distributions. These 
distributions arise from a generalized central limit theorem [29] relevant to random variables having 
an infinite variance. Specifically for n = 2, one recovers the Holtsmark distribution that describes 
the statistics of the gravitational force acting on a star in a cluster [30] and of the velocity field 



created by point vortices in turbulent flows 



3l|. 



However, for realistic non-singular flow fields, corresponding to finite values of A, we generally 
have fi > 0. Specifically, by matching the exact velocity moments to those in Eq. ©, one finds 
that for the co-oriented model (|4al) with n — 1 

c = \ N (x) 2 {kV)2 ^ 1+£+2£ 1o § ^ ( 10 ) 




FIG. 2: Velocity PDF of a tracer particle in the flow generated by different concentrations of swimmers. The 
solid curves are based on approximation (0, using the exact second and fourth moments of the velocity PDFs, 
as shown in the Supplementary Material. (A) For the co-oriented model ([4a]) with long-range hydrodynamics 
n = 1, the velocity PDFs from simulations (symbols) converge rapidly to the Gaussian distribution predicted 
by the central limit theorem (solid curves), even at low volume fractions cp ~ 0.4%. (B) In contrast, for 
the co-oriented model with n = 2, the central limit theorem convergence is very slow and the velocity 
distribution exhibits strongly no n- Gaussian features at volume fractions similar to those realized by recent 
experiments [6[. (C) The velocity PDF for the dipolar swimmer model looks very similar to that of our 
co-oriented model (B), which means the angular dependence does not play an important role for the velocity 
distribution. Simulation parameters are given in A, the sample size is 4 x 10 6 throughout. 

at leading order in I \— A/ A, whereas for n — 2 

c = KC> 1/3(Ky)V/3 ' (iia) 

For the dipolar model (l4b]h one obtains the same scaling of (c, n) with ((/?, A) as in Eq. (fTTj) but 
with a slightly different numerical prefactor, yielding in the small cutoff limit A — >> 

<M 2 > = f W 2 Kx)< (12a) 
<M 4 > = ^W 4 Kx) 5 - (12b) 

Note that Eq. (fTOl) suggests for colloidal- type flow fields with u oc r _1 the appropriate thermo- 
dynamic limit is given by TV, A — >► oc such that TV/ A 2 = const., whereas we must fix ip — const. 
if u oc r _n ,n > 2. Furthermore, Eqs. (Illbl) and fll2al) imply that /i — >> and (|v| 2 ) — > oc for a 
vanishing regularization parameter A — > 0. This illustrates that Levy-type behavior becomes more 
prominent the more "singular" the velocity field in the vicinity of the swimmer. 
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The solid curves in Fig. [2] are based on approximation 0, using the exact second and fourth 
moments of the velocity PDFs, as given in the Supplementary Material. For u oc r -2 , the Gaussian 
prediction of the central limit theorem becomes accurate only at large volume fractions (ip > 25%). 
In the dilute regime (p <C 1, the bulk of the probability comes from a Levy stable distribution before 
it crosses over to quasi-Gaussian decay, reflected by the (truncated) power-law tails in Fig. [2]B and 
C. We may thus conclude that the fluid velocity in a dilute swimmer suspension is a biophysical 
realization of truncated Levy-type random variables [^(J. 



B. Flow field autocorrelation 



The similarity of Figs. E]B and C suggests that the angular flow field structure is not important 
for the equal-time velocity distribution. By contrast, the velocity autocorrelation function CW,a(£) 
depends sensitively on the angular details, as illustrated in Fig. [3l For both our co-oriented 
model (jUD with n = 1,2 and the dipolar model (l4b]h the function Cjsf,h{t) can be determined 
analytically (see Supplementary Material). From the exact results, one finds that for n — 2 in the 
thermodynamic limit at long times t ^> r e := e/V 

C N Mt)^—^(KV) 2 (-f). (13) 

For comparison, the velocity autocorrelation function for dipolar swimmers can be approximated 

by 



1 ^ I* ^ 7 111 ' * - ' l-x A \ 

(14a) 



where £* := 41/ tt and s := tV/A < 1. The approximation f!14ap . shown as the dotted line in Fig. El 
becomes exact at long times. In the thermodynamic limit, it reduces to 



c N , A (t) * v (£) ^ { , "' ( i4b ) 

where t\ := 4\/(ttV). Note that Eq. f)14bl) predicts an asymptotic t~ 3 decay, which is considerably 
faster than the t~ x decay for the co-oriented model, cf. Eq. (fT3l) . This is due to the different 
angular structure of their respective flow fields. The excellent agreement between simulation data 
and the exact analytic curves (solid) in Fig. [3] also confirms the validity of our simulation scheme 
(see Numerical Methods in Supplementary Material). 
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FIG. 3: Due to the different flow topologies, the velocity autocorrelation function for dipolar swimmers (|4b|) 
decays faster than that for the co-oriented model (|4a|) with n — 2. Solid curves indicate the exact analytic 
solution and symbols correspond to simulation data. Dotted and dashed curves illustrate respectively the 
long time approximation and its behavior in the thermodynamic limit. Parameters are those of Fig. O 

C. Mean square displacement 



Having discussed the velocity statistics, we next analyze the tracer displacements. To that 
end, we will focus on the practically more relevant dipolar swimmer case, and include the effects 
of thermal Brownian motion (so Dq > 0). A first quantifier, that can be directly measured in 
experiments, is the mean square displacement (A[x(t)} 2 ) ~ (A[x(t)] 2 )N + 6Dot. Assuming spatial 
homogeneity and spatially decaying correlations, the velocity autocorrelation function can be used 
to obtain an upper bound for the swimmer contribution 



(A[x(t)} 2 ) 



N 



/ v /' 

JO JO 



dr (u n (t ,x{t'))u n {t,x(t))) 



< [ dr' [ dr{u N (r\ 0)u N (r,0)). 
Jo Jo 

Inserting the approximate result (Il4b|) , we find 

t < r A , 



(15) 



(A[x(t)Y) N < 6(p K 2 Ve t < 



2t 
5r A 



35r? 



t ^ ~W 



35F 



(16) 



1-^ + ^-^ t >r x . 

Equation (fT6l) implies that tracer diffusion due to the presence of dipolar swimmers is ballistic at 
short times t <C t\ and normal at large times t » t\ (see Fig. [4]). Qualitatively, the predicted 
linear growth (A[x(t)} 2 ) oc t agrees with the experimental results of Leptos et al. p. Generally, 
the asymptotic diffusion constant will be of the form D ~ Dq + vtp f^Vc, where v is a numerical 
prefactor of order unity that encodes spatial correlations neglected in Eq. (fT5l) . 
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FIG. 4: Mean square displacement at different volume fractions. Solid lines are analytic upper bounds, 
6(A) + tpK 2 Ve)t > (A[cc(£)] 2 ), and symbols display the ensemble average from simulations. Parameters are 
those of Fig. [2] 

D. Evolution of the position PDF 

The spatial motion of a passive tracer in a fluctuating flow u(t, r) is described by the position 
PDF P(t,r) = (S(r — x(t))). For Gaussian fields, uniquely defined by the two-point correlation 
function (uj{t, r)uj{t r , r 7 )), it is possible to characterize P(£, r) analytically for some classes of tra- 



jectories x(t) 



321 ]. However, our analysis above indicates that the statistics in an active suspension 



are neither 5-correlated nor Gaussian, exhibiting features of Levy processes. 



Generally, the hierarchy of correlations in Levy random fields is poorly understood [33 
therefore unclear how to adapt successful models of random advection by a Gaussian field 



. It is 



m or, 



more generally, extend the understanding of colored Gaussian noise [32] to colored Levy processes. 
These theoretical challenges make it very difficult, if not impossible, to construct an effective 
diffusion model that bridges the dynamics of P(t 5 r) on all of the time scales. Partial theoretical 
insight can be gained, however, by considering the asymptotic short and long time behavior. 

At short times t C r^, the position PDF combines ballistic transport from constant swimmer 
advection and diffusive spreading from thermal Brownian effects. For experimentally-relevant 
parameters [6], normal diffusion is much stronger than the advection and, at these times, the 
dynamics of P(t, r) are captured by the normal diffusion equation. If Brownian motion is neglected, 
we have P(t, r) = £ _3 </W,a(V /£) with 4>n,a(v) as the tempered Levy velocity PDF. The resulting 
"ballistic" Levy distribution is in good agreement with simulation data for Dq = at short times, 
see inset of Fig. [5] A. 
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For long times t ^> r\, after the correlations of the velocity field have vanished (typically after 
several seconds), we may interpret a tracer diffusing in an active suspension as a realization of 
an uncorrelated tempered Levy process. Effectively, this corresponds to replacing Ujy(t,r) from 
Eq. (pQ) with a 5-correlated but non-Gaussian random function C(t). To characterize the statisti- 
cal properties of the swimmer-induced noise C(t), we need to specify its characteristic functional 
.F[t;fe(s)] [35! ]. Our earlier findings, that the asymptotic mean square displacement grows linearly 
in time and that the velocity field amplitudes follow a tempered Levy stable distribution, suggest 
the functional form 



Here, D a is an anomalous diffusion coefficient of dimensions m a s" and the regularization param- 
eter K has dimensions m _1 . For a = 2, £(£) reduces to Gaussian white noise. Using an approach 
similar to that of Ref. the Fokker-Planck equation corresponding to Eq. (fT71) is found to be 



where we also included the contribution from normal diffusion. In Fourier space, the solution of 
Eq. (USD reads 



Eq. (IT9l) provides a good fit to the the long time simulation data in Figs. [5] C and [6J especially 
in the asymptotic regime. It is worth emphasizing that, although the motion of the tracers at 
long times is non-Gaussian and described by a fractional diffusion equation, the asymptotic mean 
square displacement exhibits normal growth, (A[x(t)] 2 ) oc t. 

On intermediate time scales, when the velocity autocorrelations are already decaying, but still 
non-neglible due to their t~ 3 scaling (see Fig. [3]), the transient behavior of the position PDF can 
interpreted as a superposition of two distinct components: i) Quasi-ballistic tracer displacements, 
which are remnants of the short time dynamics and may dominate the tails of the tracer position 
distribution, and ii) fractional diffusive behavior due to the onset of tracer scattering by multi- 
ple swimmers. A quantitative comparison suggests that the measurements of tracer diffusion in 
Chlamydomonas suspensions by Leptos et al |6[, who focused on the range t < 2 s, are exploring 
this intermediate regime. 




(17) 




(18) 




(19) 
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FIG. 5: Radial position PDF of a tracer in a dilute suspension of dipolar swimmers at various times. Solid 
curves represent analytic forms of P(t,r) and symbols illustrate histograms determined from simulations. 
Insets B and C show the same quantities on a log-log scale. Volume fractions and symbols are those of Fig.HJ 
Parameters are those of Fig. [2]with Dq = 0.245 /im 2 s _1 . (A) Short time regime. At these times, Brownian 
motion effectively dominates constant advection for our choice of parameters. In the limit of no thermal 
noise (inset), the position PDF is the tempered Levy velocity PDF 4>n,a(v) after a rescaling with t. (B) 
Transient behavior. This period corresponds to an intermediate decay of the velocity autocorrelation. (C) 
Asymptotic long time regime. Eventually, random advection from many low Reynolds number swimmers 
becomes equivalent to a tempered Levy flight. The solution to the fractional diffusion equation ([T8]) is 
matched against simulations by fitting its coefficients D a and K, see also Fig. El 
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FIG. 6: Time evolution of the marginal position PDF at a volume fraction cp = 1.6%. Solid curves represent 
Eq. (fT9j) with fitted coefficients (shown only at long times) and symbols illustrate histograms determined 
from simulations. Parameters are those of Fig. [2] with Do = 0.245 fim 2 s _1 . At intermediate times t ~ 1 s, 
our data resemble the measurements from Ref. [6[. 

III. CONCLUSIONS 



Understanding the mixing and swimming strategies of algae and bacteria is essential for de- 
ciphering the driving factors behind evolution from unicellular to multicellular life [111 ]. Recent 
experiments on tracer diffusion in dilute suspensions of unicellular biflagellate Chlamydomonas 
reinhardtii algae [6] have shown that microorganisms are able to significantly alter the flow statis- 
tics of the surrounding fluid, which may result in anomalous (non-Gaussian) diffusive transport of 
nutrients throughout the flow. 

Here, we have developed a systematic theoretical description of anomalous tracer diffusion in 
dilute, active suspension. We demonstrated analytically and by means of GPU-based simulations 
that, depending on the distance scaling of microflows, qualitatively different flow field statistics 
can be expected. For colloidal- type flow fields that scale as r _1 (due to the presence of an external 
force), the local velocity fluctuations in the fluid are predominantly Gaussian even at very small 
volume filling fraction, as expected from the classical central limit theorem. By contrast, flow fields 
that rise as r~ 2 or faster in the vicinity of the swimmer will exhibit Levy signatures. Very recent 
measurements by Rushkin et al {3^] appear to confirm this prediction. When the statistics are 
non-Gaussian, our results show that the asymptotic convergence properties of velocity fluctuations 
in active swimmer suspensions are well-approximated by truncated Levy random variables {20! ]. 

With regard to experimental measurements, it is important to note that the tracer velocity is a 
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well-defined observable only if thermal Brownian is negligible (corresponding to the limit Dq = 0). 
Otherwise, the associated displacements over a time-interval At contain a component that scales 
as y/~Kt. This fact must be taken into account, when one attempts to reconstruct velocity distri- 
butions from discretized trajectories: If thermal Brownian motion is a relevant contribution in the 
tracer dynamics, the measured distributions will vary depending on the choice of the discretization 
interval. 

Our analysis further illustrates that, for homogeneous suspensions, the angular shape of the 
swimmer flow field is not important for the local velocity distribution, which is dominated by 
the radial flow structure. By contrast, the temporal decay of the velocity correlations sensitively 
depends on the angular topology of the individual swimmer flow fields. Specifically, our analytic 
calculations predicts that velocity autocorrelations in a dipolar swimmer suspension vanish alge- 
braically as t~ 3 . This prediction could in principle be tested experimentally by monitoring the flow 
field at a fixed point in the fluid, using a setup similar to that in Ref. 

EL 

Finally, we propose that the asymptotic tracer dynamics can be approximated by a fractional 
diffusion equation with linearly growing asymptotic mean square displacement. It would be inter- 
esting to learn whether the fractional evolution of the tracer position distributions at long times 
can be confirmed experimentally. This, however, will require observational time spans that go 
substantially beyond those considered in Ref. 

In conclusion, the correct and complete interpretation of experimental data requires an exten- 
sion of Brownian motion beyond the currently existing approaches Q, Hj- Although many 
challenging questions remain - in particular, regarding the consistent formulation of a generalized 
diffusion theory that combines Levy-type fluctuations with time correlations - we hope that the 
present work provides a first step towards a better understanding of present and future experiments. 
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These supplementary notes summarize analytical calculations for the local velocity distributions, 
the velocity autocorrelation functions and the mean square displacement of tracer particles in dilute 
suspensions of active swimmers. Section [I] introduces the model and briefly summarizes a few 
well-known facts on characteristic functions required for the subsequent parts. Section HI1 discusses 
details and steps that are essential for the derivation of analytical results for the velocity statistics. 
To illustrate the general procedure, we will consider a simplified swimmer model that is helpful for 
understanding how the distance scaling of effective swimmer flow fields is reflected in the velocity 
statistics of tracers. Section IIHI provides an analogous analysis for a more realistic (regularized) 
dipolar swimmer model. Section HVl summarizes details of the numerical simulations. 



I. SUMMARY OF MODEL ASSUMPTIONS 
A. Tracer dynamics 

We consider the 3D motion x(t) = (xi(t)) of a passive tracer particle in a fluid that contains a = 1, . . . , N active 
particles (such as algae or bacteria), which are described by phase space coordinates T := {T a } = {(X a (t),X (£))}. 
We assume that, in good approximation, the tracer particle does not affect the swimmers, so that T(t) is approximately 
independent of x(t). Neglecting Brownian motion effects (corresponding to Do = in Eq. [1] of the main text), the 
tracer motion can be described by the overdamped equation (low Reynolds number or Stokes regime) 

N 

x(t) = u N (t, x) , u N (t, x) = Y, uixlT^t)). (1) 

Here, Un denotes the velocity field generated by N swimmers and u(x\T a ) the contribution of an individual swimmer a 
to the fluid flow at position x. In Sections ITT1 and Hill below, we shall study two different models for u(x\T a ). 

B. Swimmer dynamics and statistics 

We restrict our considerations to the dilute limit of small swimmer concentrations. In this case, we may assume 
that the swimmers move approximately on straight lines at constant velocity X (t) = Vq. Then, T a (t) is uniquely 
determined by its initial condition Tq := T cr (0) = (Xq , Vq), and we have 

X%t) = Xl + tVZ- (2) 

Throughout, it is assumed that the swimmer speed is roughly the same and equal to V and that, initially, the swimmers 
are uniformly distributed in a large spherical container (radius A), i.e., the initial conditions Tq = (Xq, Vq) are i.i.d. 
random variables with joint probability density function (PDF) 

<Mro) = ll*i(r ), $ i( r o)= (47r/3)A3 4wV 2 > ( 3 ) 
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where 3(y\V) is the indicator function of the spherical container volume V = {x : \x\ < A}, 

J( V |V) == Z 1 ' V ! V ' (4) 

We are interested the velocity distribution at the center of the volume, formally defined by 

^n,a( v ) := (<S(«-ujv(0,O))), (5a) 

and, moreover, in the flow field autocorrelation function at a fixed point in the fluid far from the boundaries, 

CW,a(*, ) := (u N (t, 0)u N (0, 0)), (5b) 

where ( • ) indicates an average with respect to the initial swimmer distribution (j3]) . Generally, we are interested in 
evaluating these quantities in some suitably defined thermodynamic limit N, A — >> oo. 

C. Characteristic function 

Consider the characteristic function \$ °f some velocity PDF <fr(v), defined by 

X4>(0) : = J d dv exp(-iqv) <j>(v). (6) 
Then, by inverse Fourier transformation, we have 

^ = J^kY j exp (^ v ) X4>{0)' ( 7 ) 
Given x</>(<7), the second and fourth velocity moment can be obtained by differentation 

<M 2 > = -A qX0 (O) , (| V | 4 ) = A^(0). (8) 

In particular, for 3D spherically symmetric velocity PDF <j>(y) = (j)(v) with v = \v\, we have x</>(q) — X</>((Z)> wnere 
q = \q\. In this case, the differential operators simplify to 



A q = ^ + -d q j , = ^ + -d q j = d q + -d*. (9a) 

Below, we shall show that the velocity PDF of a tracer in the presence of active swimmers can be approximated by 
a tempered Levy stable distribution with characteristic function 

XM = exp{-[(c q 2 + n 2 ) a/2 ~ » a }}, (10) 
which yields for the second and fourth velocity moment 

(|i;| 2 ) = 3ac//*- 2 , (11a) 
(\v\ 4 ) = 15ac 2 /i a_4 [2 + a(/i a -l)]. (lib) 

The tempered Levy distribution in Eq. (fTO]) reduces to either the Gaussian prediction from the central limit theorem 
as \± becomes large (many swimmers) or the Levy stable distribution from the generalized central limit theorem [1] 
as \i — >• (unregularized swimmers). 

II. FIRST EXAMPLE: REGULARIZED CO-ORIENTED MODEL 

We first consider the model interaction 

u(x\T5) = (kV%) ! — -4r r"> ( 12a ) 

where e can be interpreted as the swimmer radius, < n < 1 is a dimensionless coupling constant, A a regularization 
parameter, Vq the swimmer velocity (assumed to be constant), and 

r a (t) :=X a (t)-x, (12b) 

is the vector connecting the swimmer position at time t, X a (t) = Xq + tV^ , with the space point x. 
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A. Velocity PDF of the tracer particle 

Inserting the Fourier representation of the ^-function, we can rewrite ([5a]) as 

«-$>(o|rg) 



(13) 



Using spherical variables, Xq = RQ and Vq = VCl for both the initial swimmer positions and velocities, we find 



(47T/3)- 



A 3 



^3 J tfqeMidv) j^AR R 2 J d 2 Q exp(-i kV J^^Q^ 



N 



(14) 



Introducing a rescaled radial position variable y = ii/A, and performing the angular integral over d 2 ^, we obtain 



</>n,a( v ) 



(2tt; 



^3 y 4 d 3 gexp(zqrv) if^(g), 



where 



is the characteristic function of < 



jtf( g ) = fa / dy , 2Sin[A " (2/)l9l] 

L ^0 

with 



My)\i\ 



N 



(e/A) n 

A n (y) := 1 L-L- kV. 
It is useful to rewrite in the equivalent form 

r r 1 r 1 i N 

K%{q) = 3/ df/ dyy 2 cos(A n \q\0 =: [K n (q)] N , 
. Jo Jo 

where K n (q) is the characteristic function for the one-swimmer case. 



(15a) 



(15b) 



(15c) 



(16) 



1. Exact second and fourth velocity moments 



We note that, for m € N, 



K n (0) = 1, 
d 2m K n (Q) = (-iy 
dj, 2m -^K n (0) = 0, 



2m + 1 Jo 



dy y 2 A 2m , 



\\m -df m -^K n {q) = (-1)™-JL- (\yy 2 A 2m = d 2m K n (0). 
q^o q H zm + 1 Jo 



From these expressions, one finds 



and, therefore, 



d q K?(0) = 0, 

%K?{0) = N[d 2 K n (0)}, 

d 3 q K»(0) = 0, 

%KZ(0) = 3(AT - l)N[d 2 K n (0)} 2 + N[%K n (0)], 



AK»(0) = 3N[d 2 K n (0)\, 

A 2 K%(0) = lhN(N -l)[d 2 K n (0)} 2 + 5Nd 4 q K n (0). 



(17a) 
(17b) 
(17c) 
(17d) 



(18a) 
(18b) 
(18c) 
(18d) 

(19a) 
(19b) 
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Equations (fT9j) hold for any spherically symmetric distribution in 3D. 

From the above expressions, one finds the following exact formulas for the velocity moments in this model: 

<M 2 ) = 3N (\yy 2 [A n (y)} 2 , 



4 5(N-1) 

where, using the abbreviation t := A/A, 

- 1 (kV) 2 /e\2n r r 



\v\ 2 )^ + ZN f&yy 2 [A n (y)}\ 
Jo 



J dyy 2 [A n (y)} 2 = { —L (_) |__ + _(„ _ 3)HypGeom2Fl 



3(n- 1) 



(2n - 3) r + (5n - 6) £ 2n + 



i 3 3+n _ ' 



2n , 9 - 18n + lln 2 



?3n 



(1+^ 

(n - l)(n - 3)(2n - 3)HypGeom2Fl 



3(n- 1) 
3 3 + n 



}• 



For example, for the (n = 2) model, these expressions simplify to 







dy y 2 [A n (y)Y 



dy y 2 [A n (y)} 4 



1 /e\ 4 
2 



arctan 



l6~(a) (a) {(? 



arctan 



1 



l + £ 2 
1 



V (l + £ 2 f 3(1 + £2)3 



(20a) 
(20b) 

(20c) 
(20d) 

(21a) 
(21b) 



The first integral, and hence (|i;| 2 )0, diverges as 1 / A if one lets the cut-off A — >• (which is equivalent to £ —> 0). More 
precisely, in this limit, Eqs. ([20]) reduce in leading order of (A/ A) to 



<M 2 >« 
(\v\ 4 ) 4 



3tt 

T 

3tt 
32 



(«v)'„(i) 
(«y) 4 ^ ( 



ex5 
A/ 



(22a) 
(22b) 



where 9? := A?"(e/A) 3 is the volume filling fraction. 



2. Approximation by a tempered Levy distribution 
We would like to approximate the exact characteristic function from Eq. (|16p by the tempered Levy law 

X*{q) = exp{-[(c q 2 + M 2 )"/ 2 - (23) 
which exhibits quasi-Gaussian behavior for small q, corresponding to large velocities, 

XM ^ exp (-afi a - 2 c . (24) 

To motivate the ansatz (|23]h we write (|15b|) as 

K^(q)=exp{-N\n[K n (q)}}, (25) 

and consider the limit A = 0. In this case, the double integral for the one-swimmer characteristic function K n (q) 
from Eq. ([T6]) can be calculated exactly for n = 1,2,3 in terms of trigonometric, hypergeometric, and sine integral 
functions. By expanding the resulting expressions \n[K n {q)) for a large volume A ^> e, one obtains for n = 1 a 
Gaussian limiting distribution 

K»(q) ~ exp (-7\ q 2 ) , T,= l -N (^f (kV) 2 . (26a) 
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By contrast, for n > 2 the limiting distribution is found to be of the Holtsmark type, i.e., 

K»(q) ~ exp (-T n q 3 / n ) , T n =t n N («) 8 ( K Vf' n 

87r/5 and ts = 7r/4. For comparison, if we let /i 



where t n is a constant of order unity; specifically, £2 
Eq. ([23]) , we obtain 



X4>(q) = exp (-c a/2 g a ) , 
Thus, by comparing with ([26]) and (|26cp . we can identify 

n = 1 = 
n > 2 => a = 



a = 2, 
3 



n 



(26b) 
^ in 

(26c) 

(27a) 
(27b) 



Before determining the remaining parameters (c, /i), it is worthwhile to note that the effective temperature scales 
differently with volume and swimmer number for n = 1 and n > 2, respectively: In the case of a colloidal- type 
velocity field with n = 1, the effective temperature Ti is proportional to the area filling fraction 7V(e/A) 2 , see 
Eq. (|26ap . whereas for swimmer- type flow fields with n > 2 the effective temperature T n scales with the volume filling 
fraction 



V := iV(e/A) 3 



(28) 



see Eq. ([26b)) . This suggest that, for n = 1, the appropriate thermodynamic limit corresponds to A/", A — ^ oo while 
keeping the ratio N/ A 2 fixed, whereas for n > 2 one should let A/", A — >• oo such that Af/A 3 remains constant. 

It remains to discuss how to identify the parameters (c, /i), which we shall do next for the cases n = 1,2,3. For 
n = 1 the procedure is rather straightforward; in the case of n = 2,3, we are going to determine (c, /i) by matching 
the second and fourth velocity moments of the tempered Levy ansatz ([23|) . which were given in Eq. ([TT]h to the exact 
moments ([20]) . 

n=l.- In this case, according to Eq. ([27|). we have a = 2, and the ansatz ([23]) reduces to the Gaussian 



X0(g) =exp(-cg 2 ). 
By comparing with ([26ap . we see the c = Ti, i.e., 

2 



= ^(l) 2 ( * y)2 [i + ' + ^M*)]- 



(29a) 



(29b) 



where we now also included the leading order corrections in £ = A/ A. 

n = 2.- In this case, we have a = 3/2 and the velocity moments ([TT]) of the tempered Levy ansatz ([23]) take the 
form 



<M 2 > 
<M 4 > 

Solving these equations for (c, /i) yields 



9 
2 
45 



2 C/i 



-1/2 



cV 5 / 2 



2 /5 



9 V3 



1/3 



<M 2 > 



2\5/3 



(3(H4)-5(b|2)2) 



1/3 ' 



/i 



2/3 



<|«| 2 ) 4 /3 



(3(b|4)-5(H2)2) 



2\2\ 2 / 3 ' 



(30a) 
(30b) 

(30c) 



Here, we can inserting for (\v\ 2 ) and (|v| 4 ) the exact expressions ([20]) ; this gives the fit curves shown in Fig. 2 of the 
main paper. Furthermore, by expanding the resulting formula for large volume A ^> e, A, we find 



-Or) 



2/3 



7' - 



2/3 



(30d) 
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n = 3.- In this case, we have a = 1 and the general expression for the velocity moments of the tempered Levy 
law from Eq. ([TT]) reduce to 

3c/i _1 , (31a) 



v 



|2\ 



(|i;| 4 ) = 15c 2 /i" 3 (/i + 1) • (31b) 

Solving these equations for (c, /i) yields 

c= 5»! 5 (l^l 2 ) 2 ( 3ic) 

Inserting for (|v| 2 ) and (|v| 4 ) the exact expressions (|20|l and expanding for A 3> e, A, we obtain 



c =-(/d/)V, ^ = 5 (jJ ^ ( 31d ) 

Note that for a vanishing A — » 0, the parameter /i goes to zero in Eqs. (|30d|) and (|31d|) . which implies that the 
second moment (\v\ 2 ) diverges. This also illustrates why for n > 2 - or, more precisely, for n > 3/2 if one allows for 
non-integer exponents n - there is no convergence to a Gaussian distribution in the limit A —> . 



B. Velocity autocorrelation function 

We are interested in the fluid's velocity autocorrelation function (|5b[) for the power-law model ([T2]) . We start from 
Eq. (|5b|) . which can be written as 

c n ,aV) = ^(^(^|r a (t))^(o|rg)). (32) 

Assuming, as before, that the initial swimmer position and velocities are distributed according to (j3j), we can simplify 



CsAt) = N(kV) { (pcg + tvg|n + Am) (jx^p+A™)/ =: N{KV) ^ (33) 

Using spherical velocity and position variables, Xq = XCl and Vq = Vti, and inserting the one-swimmer PDF from 
Eq. (|3]), we find 



3 r <m i i 

4tt J ~4tt (y 2 + s 2 + 2syz) n / 2 + ^ n ?/ n + ^ n ' 



By virtue of the identity 



Eq. (|36|) can expressed as 



where 

Tic n \ • — I rl 

(?/ 2 + s 2 + 2s?/z) n / 2 + ^ n ' 



(34) 



Introducing rescaled variables 

2/ := X/A , 5 := tV/A , £ := A/A , * := fill (35a) 

and noting that 

|xn + tyn| = a [ y 2 + s 2 + 2^)] 1/2 , (35b) 

we can rewrite (l34l) in the form 



M=(j) 2n J^l'dyy*^ I ^ , , o ; N „„ , „, (36) 



JdnjdCl f(z) = (47r)(27r) y d* /(«), (37) 

CnW = 2 (a) L dy rT¥ Jn(s ' y) ' (38a) 



./„(.<•//) := | ( <U: — , „' , . .. , • (38b) 



We next provide explicit results for n = 1 and n = 2. 
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n = 1. 



In this case, we find 



M s >y) 



sy 



— log 

sy 



£+\s-y\ 



y + l\ 



(39a) 



The remaining integral integration in Eq. (j38a|) can be easily computed numerically. However, the correlation 
function c\{t) can also be expressed analytically in terms of the polylogarithm Li q (z), defined by 



k=l 



(39b) 



and analytic continuation for \z\ > 1. Li g (z) is real valued for real z < 1 and possesses the integral representation 



Li g (z) 



r(g) 7 



dfc 



(39c) 



Considering a sufficiently large volume such that < s < 1, we obtain 



ci(t) 



'I (x) 2 ^{ s2 + ( 2 ^- 2 ) 5 + 2 ^( 1 + ^) arctann (Y^ 



t 2 loe 



^ 2 



Li 2 



4 + log 
— Li 2 



1 + * 



log 



(l + £) 4 (l + £) 2 



Li 2 



(5 + £) 4 (S + £) 4 



(39d) 



Li2 



-2£ 



Li2 



'*)]}■ 



with 3? denoting the real part. In particular, in the thermodynamic limit TV, A — >• oo such that cp = N(e/ A) 2 = const., 
we find that the autocorrelation function Cjv,a(£) = N(nV) 2 c n (t) becomes constant 

C^\t) = 3 V ( K V) 2 . (39e) 

This situation, however, is unrealistic for real swimmers, which typically generate flow fields that decay with n > 2. 
n = 2.- In this case, we find 



2s?/ 



1 



4s?/ 



(s-iy) 2 +^ 2 

Considering again < s < 1, the correlation function c 2 (£) may be written in terms of the Dilogarithm Li 2 (z) as 



(40a) 



c 2 (t) 



Li 2 



log(s) lo; 

s-l-ii 



1 



4s 



Li 2 



(s-l) 2 +£ 2 
's + l-il' 



Li 2 



s + 1 + ^ 
s - 1 - it 



■Li s 



s + 2i£, 

V s-2^ y j ' 



s / v ^ y v « - J ~* V s - 2U 1 { ' (40b) 

with 5ft denoting the real part. In the thermodynamic limit N, A — >• oo such that cp = A^e/A) 3 = const., we find at 
large times t ^> r e := e/V for the full autocorrelation function 

3tt 2 



(40c) 



III. SECOND EXAMPLE: REGULARIZED DIPOLAR SWIMMER MODEL 



Let us now consider a more realistic dipolar swimmer flow field model, denned by 



v(x\r° ) = (kV) 



(n r ) 



i 



|iT| 2 + A 2 

where k > (k < 0) correspond to an extensile (contractile) swimmer, and 

R(t) 



R , 



A > 0, 



R a {t) := X a {t) -x, R (t):= 



\R (t)\ 



(41a) 



(41b) 



As before, we assume, for simplicity, that a swimmer's orientation does not change over time. The dipolar swimmer 
model (jHJ) exhibits the same distance scaling as the regularized power-law model from Eq. ([T2]) with n = 2, but the 
directional dependence is different. As a consequence, as we shall see below, the velocity PDFs of the two models are 
very similar but the correlation functions show qualitatively different behavior. 
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A. Tracer velocity PDF 

1. Characteristic function 

Using spherical variables, Xq = Xtt and Vq = Vti for both the initial swimmer positions and velocities, the 
characteristic function of the velocity PDF can be written as 



K?{q) = [K s (q)] N 
where the one-swimmer characteristic function is given by 



Ks{q) = (47r ff 1 _L [ A dRR 2 I d 2 n I d 2 n cx P 



A 3 47T J 

Hence, with y := R/A, I := A/A and z := ftfl 



j-i(/c 



V) 



X 2 + \ 2 



(qSl) 



K S (q) 



I [ dy y> J d»n j\ d* exp {- *(W) (£)° {-^) (Sz 2 - l) (,11) 



Choosing q = (0, 0, \q\) = (0, 0, q), with no loss of generality, we find 

K S (q) = \f Ayy2 j\ dz j Q d ^ cos (^)' 



where 



D(y,z) := (kV) Q)' 



(3z 2 - 1) . 



„y 2 + P, 

This result can be used to evaluate analytically the moments of the tracer velocity PDF. 



(42) 



(43) 



(44) 



(45a) 



(45b) 



2. Velocity moments 

From Eqs. (jIS"|) and (fT§|l . we find the following exact results for the second and fourth velocity moments: 

arctan (] 



(M 2) = N{kV) 2 Q) 4 j6 

<M 4 > = N(KVf (1) 



£ 



i+^j / 



3 

35 



3arctan(l/f) 3 + S£ 2 - 3£ 4 



£ 5 



where 



A/ A is the rescaled regularization cut-off. In the small cutoff limit 

3tt 



57V- 1 
+ 3 N 

0, we find 



(\v 



1 2 \ 2 



(46a) 
(46b) 

(47a) 
(47b) 



These expressions are quite similar to those obtained for the power-law model with n = 2, see Eq. (|22]h The moments 
(|47j) can be used to determine the parameters of the corresponding tempered Levy velocity distribution ([23]) by means 
of Eqs. (|30dj) . However, as we shall see below, the two models give rise to very different velocity correlations. 



B. Velocity autocorrelation function 



We are interested in the velocity autocorrelation function (|5b|) of the fluid near the center of the volume, which can 
be written as 



C s (t)=N(KV) 2 



R(t)\ 2 + \ 2 J V|i^(0)| 2 + A 2 



QR(t) 



ClR(0) - 1 \ R(t)R(0) 



(48a) 
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For linear swimmer motions, we further have 

R(t):=x + tv , x = xn, v = vn 



(48b) 



with X and V denoting the initial swimmer position and velocity; ft and Cl are the corresponding unit vectors. 
Defining the dimensionless velocity autocorrelation function cs(t) by Cs(t) = N(kV) 2 cs(t), we obtain 



cs(t) 



\X + tVn\ 2 + \ 2 J Vl^l 2 + A 2 



/ xnn + tv" 
y\x + tvn\ 



x + tvnn \ 
\x + tvn\ / ' 



(49) 



In the limit t = 0, one recovers the second velocity moment, cs(0) = (\v 2 \) / [N (nV) 2 ]. 
For t > 0, using the notation from Eqs. ([34]) and (|35|h we have 



cs(t) = e 4 



1 



3^-1 



3(Xz + tV) 2 



X 2 + t 2 V 2 + 2tVXz 



1 



K X 2 + \ 2 J \X 2 + t 2 V 2 + 2tVXz + A 2 
Substituting y := X/A, t = A/A, s := tV/A and using the identity ([57)1 . this can be written as 

3 /e\ 4 



VX 2 + t 2 V 2 + 2tVXz 



c s (t) = - (1) AJ(t) , AJ(t) := Ji(t) - / 2 (i), 



where 



dz 



3z 2 -l 



£ 2 + y 2 J \£ 2 + y 2 + s 2 + 2syz 



y + sz 



h{t) = dy dz 



y 



3z 2 -l 



£ 2 + y 2 ) \£ 2 + y 2 -f s 2 + 2syz 



(y 2 + s 2 +2syzy/ 2 _ 
y + sz 



3(yz + s) 2 



_{y 2 +s 2 + 2syz) 2 l 2 _ 



(y 2 + s 2 + 2s^)V2 



The integration over z can be performed analytically, and the remaining ^/-integrals can be written as 



h(t)= dyh&sj), h(t)= dyj 2 (y;s,£) 



where 



3 (3s 4 + 2s V + 3y 



32s V^ 2 (y 2 +£ 2 ) 
3£ 2 (s 2 - 9y 2 - h£ 2 ) 



) 

J- [(s - yf(s + y) 2 - (s + y) 3 (s - y)\s - y\ 



[{s + y) + \s-y[ 



80s 4 y (y 2 +£ 2 ) 
3£ 2 [-32s 4 + 18y 4 + 40y 2 £ 2 + 15£ 4 + 2s 2 (19?/ 2 + U 2 )] 

160s V (y 2 +£ 2 ) 
201s 4 + 185s V - 174y 4 



[(s + y) - \s-y[ 



280s 5 ?/ (y 2 +£ 2 ) 
201s 4 + 596s V - 381y 4 
560sV (y 2 +£ 2 ) 



[(s + y) 2 + (s - y)\s - y\] + 
[(s + ^) 2 -(s-^)|s-^|] + 



3 (s 2 - ?/ 2 - £ 2 ) 2 3s 6 - s 4 (?/ 2 - 9£ 2 ) + 8s 2 £ 4 + (?/ 2 + £ 2 ) 2 (s 2 - 3?/ 2 + 3£ 2 ) 



32s V^ 3 (?/ 2 +^ 2 ) 
s + ?/ 



arctan 



£ 



arctan 



I s -2/1 



J2(y;s,£) 



(y 2 - 9s 2 - 5^ 2 



20s 2 ?/(7/ 2 +^ 2 ) 

33s 4 - 177/ 4 + 10y 2 £ 2 + 15£ 4 + 8s 2 (?/ 2 + 5£ 2 ) 
40s V (y 2 +£ 2 ) 

(s 2 - t/ 2 + £ 2 ) [3s 4 + 3 (?/ 2 + £ 2 f + 2s 2 (y 2 + 3^ 2 ) 



[(s + 2/) - |s-2/[ 



8s 3 y 2 £(y 2 +£ 2 ) 



arctan 



s + y 



arctan 



\ s -y\ 



(50) 

(51a) 

(51b) 
(51c) 

(52a) 



(52b) 



(52c) 
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To obtain the exact correlation function, the remaining one-dimensional ^/-integrals (|52a|) can be computed numeri- 
cally; for special limit cases, however, one can expand the integrands .71/2(2/; s,i) and evaluate the resulting integrals 
analytically. 

Short-time expansion t —> 0. Expanding the integrands .71/2(2/; s,£) at short times s <C we find 



C s (t) ~ v ( 



176s 5 



66s 3 



+ ^ 2 (s 2 + £ 2 ) 3 (s 2 + Pf 



(489 + 1088£ 2 + 423£ 4 ) 
<? 2 (l + £ 2 ) 3 



489 

T 2 " 



237 \ 3 (163s 2 - 84£ 2 ) [arctan(€) - arctan (£/s)] 



£ 2 



}■ 



(53) 



where £ A/A and s := £V/A. 

Large-time (small £) approximation and thermodynamic limit. To obtain an analytically tractable approximation 
of the autocorrelation function that can be used to determine the thermodynamic limit, we note that, at large times, 
we can approximate £ ~ in the denominators of the integrals Eqs. (|51b|) and (|51b|) . With this simplification, the 
^-integral can be computed more easily, yielding for < s < 1 



16 (7s 2 y - Qy 3 
35s 5 



9(s - y) + 



8 (-9s 2 + 7y 2 
3%4 



Q(y-s), 



(54) 



where 6(x) is the unit step function, defined by Q(x) := 0,x < and Q(x) := l,x > 0. In principle, the remaining y 
integral can readily be evaluated to obtain the long-time behavior of Cs(t). However, the resulting expression diverges 
at short times, since we let the cut-off £ — » 0. To avoid this divergence and mimic the effect of the short-distance 
cut-off £, we may replace the lower integral boundary in Eq. (p3|) by a regularization parameter £* [3[, i.e., we compute 



Cs(t)^ {K V?(^j\y 



16 (7s 2 y - 6y 3 
35s 5 



e(s-y) + 



8 (-9s 2 + 7y 2 
35? 



G(y - s) 



(55) 



which gives 



3s z 
7 



7s 5 



3s z 
711 

7 



£* < s < 1. 



(56a) 



The regularization parameter £* can be determined from the condition Cs(0) = (|v| 2 ), by using the exact result ()46a|) 
for the second moment, yielding 



£* - -i 

7T 



(56b) 



The result (|56|) becomes exact at sufficiently large times s 1; it does, however, also provide a useful approximate 
description at intermediate and small times. In particular, the second line in Eq. (J56a| implies that, for a finite 
system, the autocorrelation function becomes negative after a certain time to, which can be estimated as 



to 



£<i f 4 A 

7T A 



2/3 



(57) 



Physically, this due to the dipolar flow field structure: If a dipolar swimmer passes a fixed point in the fluid, the flow 
at this point will reverse its sign (direction) after certain time. 

However, as evident from Eq. ([57]) . the negative correlation region vanishes for A — >> 00, since the zero to of Cs(t) 
moves to 00 in this limit. More precisely, by taking the thermodynamic limit of Eq. (|56|) . we find that 



Cs(t) 



t<T X :-- 




1 A 

7T V 



(58) 



Thus, the velocity field autocorrelation function in a dipolar swimmer suspension decays asymptotically as t 3 , which 
is different from the t~ x -decay found earlier for the power-law model with n = 2, see Eq. (|40c[) . 
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C. Spatial mean square displacement of a tracer particle 

The approximate result for the velocity autocorrelation can be used to obtain an upper bound for the swimmer 
contribution to the mean square displacement (MSD) of the tracer particles. Considering an initial tracer position 
x(0) = and noting that (un(t' ', 0)uat(t, 0)) = C{\r' — r|), we find 

(lx(t)} 2 ) = f dr' f dr (u n (t',x(t'))u n (t,x(t))} 
Jo Jo 

rt r t 



< f dr' f dr (u n (t',Q)u n (t,Q)) 
Jo Jo 



t pt 

dr' / drC{\r' -t\). (59) 
o Jo 

The second line reflects, roughly speaking, the assumption that correlations are spatially homogeneous and decaying 
with distance, i.e., 

{u n (t', x')u n (t, x)) < (u n (t \x)u n (t, x)) = (u n (t\ 0)u n (t, 0)) (60) 
Changing integration variables, r \-> 6 := r' — r, and, subsequently, r' ^ v := r' — t, we may rewrite Eq. (f59j) as 

([x(t)] 2 ) ~ j\£ +t d6C(\9\). (61) 

Inserting the approximate result ([58]) for the fluid autocorrelation function in the thermodynamic limit, we find 

([x(t)f) c 6^ 2 Ve J^"^ 2 ; 2 t4 t<Tx:= ^> (62) 

I 1 t ^ 5t 2 350 ' 1 — Tx ' 

Thus, tracer diffusion in a dilute suspension of dipolar swimmer is ballistic at short times t < ta and normal at large 
times t^> t\. 

IV. NUMERICAL METHODS 

In our computer simulations, we directly integrate the Langevin equations 



x(t) = u N (t,x(t)) + y/2Do£(t), (63a) 

N 

u N (t, x) = ^l r " (*)) > r " W = ^ a W), (63b) 

cr=l 

X ff (t) = X°+tV° , , a = l,...,N, (63c) 

where £(t) is Gaussian white noise, to obtain the velocity distribution at a given point in the fluid and velocity 
autocorrelaton functions (for Dq = 0), and the tracer position PDF (Dq > 0). The initial positions and velocities 
r°"(0) = (Xq, Vq) of the swimmers were sampled from the distribution ((3]). 

Particle deletion & insertion.- Using the Euler method, we simulate an "ideal gas" of active particles (swimmers), 
which move according to Eq. (j63cp through sphere of radius A. This sphere is always relative to the passive tracer, 
whose position evolves according to Eq. (j63aj). If an active particle leaves the sphere, we immediately delete it. 
To restore detailed balance, we continually insert new active particles at the boundary of the sphere. The number 
of insertions per time step is drawn from a Poisson distribution f(j) = v 3 e~ v / j!, where v is the mean number of 
insertions during At. We obtain equilibrium by setting v equal to the mean number of deletions during At, which 
may be estimated from the kinetic theory of gases [2| as v = ^j^-At. For each insertion, it is necessary to bias 

the orientation of an active particle so that its probability distribution satisfies p(tl \R ) oc Cl R normalized over 
the solid angle of a hemisphere with inward surface normal R . We achieve that by uniformly choosing a position 
R a at a distance A relative to the tracer, then choosing an orientation from p(£l \R ). Generally, our procedure 
achieves numerical accuracy if A is large, and ensures that a suspension of ballistic particles remains homogeneous 
and isotropic with mean population N. A comparison with the exact results for time-dependent velocity-correlations 
verifies the validity of this approach. 
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GPU implementation.- Resolving the tail of a probability distribution can be a computationally expensive task, 
even for stochastic processes that are relatively simple. In our numerical calculations, further difficulties arise from 
having to create and maintain an active suspension unique to each tracer. This process would not be possible in a 
reasonable amount of time on a traditional computer. We therefore implemented parallelized simulations on a graphics 
processing unit (GPU) using NVIDIA's Compute Unified Device Architecture (CUDA). Compared to a single-core 
CPU, GPU code yields substantial speed-ups (up to a factor of a few hundreds). However, our longest simulation of 
4194304 trajectories still took 14 days on a GPU. 
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